// --*- C++ -*------x-----------------------------------------------------------
//
//
// Class:           SWAlign
//
// Author:          Enrico Negri, Alessandra Mantovan, Silvio Tosatto
//
// Project name:    BiocompUP Align
//
// Date:            12/2007
//
// Description:     Implement Smith-Waterman local alignment.
//
// -----------------x-----------------------------------------------------------

#include <SWAlign.h>

namespace Biopool
{

// CONSTRUCTORS:

SWAlign::SWAlign(AlignmentData *ad, GapFunction *gf, ScoringScheme *ss)
	: Align(ad, gf, ss)
{
	pCalculateMatrix(true);
}


SWAlign::SWAlign(AlignmentData *ad, GapFunction *gf, ScoringScheme *ss,
	const vector<unsigned int> &v1, const vector<unsigned int> &v2)
		: Align(ad, gf, ss)
{
	pCalculateMatrix(v1, v2, true);
}


SWAlign::SWAlign(const SWAlign &orig) : Align(orig)
{ }


SWAlign::~SWAlign()
{ }


// OPERATORS:

SWAlign&
SWAlign::operator = (const SWAlign &orig)
{
	if (&orig != this)
		copy(orig);
	return *this;
}


// PREDICATES:

void
SWAlign::getMultiMatch()
{
	Traceback tb = B0;
	int i = tb.i;
	int j = tb.j;
	pModifyMatrix(i, j);

	tb = next(tb);
	while (((tb.i >= 0) || (tb.j >= 0)) && ((i != tb.i) || (j != tb.j)))
	{
		ad->calculateMatch(i, tb.i, j, tb.j);
		i = tb.i;
		j = tb.j;
		pModifyMatrix(i, j);
		tb = next(tb);
	}

	pCalculateMatrix(false); // recalculate F, B and B0 from modified data
	ad->getMatch();
}


// MODIFIERS:

void
SWAlign::copy(const SWAlign &orig)
{
	Align::copy(orig);
}


SWAlign*
SWAlign::newCopy()
{
	SWAlign *tmp = new SWAlign(*this);
	return tmp;
}


// HELPERS:

void
SWAlign::pCalculateMatrix(bool update)
{
	int maxi = n;
	int maxj = m;
	double maxval = INT_MIN;

	for (int i = 1; i <= static_cast<int>(n); i++)
		for (int j = 1; j <= static_cast<int>(m); j++)
		{
			double s = ss->scoring(i, j);
			double extI=0, extJ=0;

			if ((i != 1) && (j != 1))
			{
				if (B[i-1][j].j == j)
					extI = F[i-1][j] - gf->getExtensionPenalty(j);
				else
					if (B[i-1][j].j == (j - 1))
						extI = F[i-1][j] - gf->getOpenPenalty(j);
			}
			else
				extI = F[i-1][j] - gf->getOpenPenalty(j);

			if ((i != 1) && (j != 1))
			{
				if (B[i][j-1].i == i)
					extJ = F[i][j-1] - gf->getExtensionPenalty(j);
				else
					if (B[i][j-1].i == (i - 1))
						extJ = F[i][j-1] - gf->getOpenPenalty(j);
			}
			else
				extJ = F[i][j-1] - gf->getOpenPenalty(j);

			double z = F[i-1][j-1] + s;
			double val = max(max(max(z, extI), extJ), 0.00);

			if (update)
				F[i][j] = val;

			if (EQUALS(val, 0))
				B[i][j] = Traceback::getInvalidTraceback();
			else
				if (val > 0)
				{
					if (EQUALS(val, z))
						B[i][j] = Traceback(i-1, j-1);
					else
						if (EQUALS(val, extJ))
							B[i][j] = Traceback(i, j-1);
						else
							if (EQUALS(val, extI))
								B[i][j] = Traceback(i-1, j);
							else
								ERROR("Error in SWAlign: SW 1", exception);
				}
				else
					ERROR("Error in SWAlign: SW 2", exception);

			if (val > maxval)
			{
				maxval = val;
				maxi = i;
				maxj = j;
			}

			B0 = Traceback(maxi, maxj);
		}
}


// SSEA variant
void
SWAlign::pCalculateMatrix(const vector<unsigned int> &v1,
	const vector<unsigned int> &v2, bool update)
{
	unsigned int minL = 0;
	

	int maxi = n;
	int maxj = m;
	double maxval = INT_MIN;

	for (int i = 1; i <= static_cast<int>(n); i++)
		for (int j = 1; j <= static_cast<int>(m); j++)
		{
			// start SSEA variant code
			if (v1[i-1] < v2[j-1])
				minL = v1[i-1];
			else
				minL = v2[j-1];

			double s = ss->scoring(i, j) * minL;
			// end SSEA variant

			double extI=0, extJ=0;

			if ((i != 1) && (j != 1))
			{
				if (B[i-1][j].j == j)
					extI = F[i-1][j] - gf->getExtensionPenalty(j);
				else
					if (B[i-1][j].j == (j - 1))
						extI = F[i-1][j] - gf->getOpenPenalty(j);
			}
			else
				extI = F[i-1][j] - gf->getOpenPenalty(j);

			if ((i != 1) && (j != 1))
			{
				if (B[i][j-1].i == i)
					extJ = F[i][j-1] - gf->getExtensionPenalty(j);
				else
					if (B[i][j-1].i == (i - 1))
						extJ = F[i][j-1] - gf->getOpenPenalty(j);
			}
			else
				extJ = F[i][j-1] - gf->getOpenPenalty(j);

			double z = F[i-1][j-1] + s;
			double val = max(max(max(z, extI), extJ), 0.00);

			if (update)
				F[i][j] = val;

			if (EQUALS(val, 0))
				B[i][j] = Traceback::getInvalidTraceback();
			else
				if (val > 0)
				{
					if (EQUALS(val,z))
						B[i][j] = Traceback(i-1, j-1);
					else
						if (EQUALS(val, extJ))
							B[i][j] = Traceback(i, j-1);
						else
							if (EQUALS(val, extI))
								B[i][j] = Traceback(i-1, j);
							else
								ERROR("Error in SWAlign: SW 1", exception);
				}
				else
					ERROR("Error in SWAlign: SW 2", exception);

			if (val > maxval)
			{
				maxval = val;
				maxi = i;
				maxj = j;
			}

			B0 = Traceback(maxi, maxj);
		}
}

} // namespace
